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Formation of stationary 3D wave patterns generated by a small point-like impurity moving through 
a Bose-Einstein condensate with supersonic velocity is studied. Asymptotic formulae for a stationary 
far-field density distribution are obtained. Comparison with three-dimensional numerical simula- 
tions demonstrates that these formulae are accurate enough already at distances from the obstacle 
equal to a few wavelengths. 

PACS numbers: 03.75.Kk 



I. INTRODUCTION 

As is well known, superfluidity means that slow enough flow of a fluid is not accompanied by heat production or 
generation of excitations of any kind. As a result, movement of a fluid is free of dissipation. In a similar way, motion 
of impurity through a superfluid goes on without any friction for small enough values of its velocity. The threshold 
velocity above which superfluidity is lost is determined by various physical mechanisms depending on the nature of 
the fluid and geometry of the process. For example, in the original Landau theory [H,[2| of superfluidity it breaks down 
when the generation of rotons becomes possible which leads to a famous Landau criterion for superfluidity. However, 
Landau's estimate for this mechanism of dissipation gives too large threshold velocity for Hell, and this disagreement 
with the theory was explained by Feynman [3| by the possibility of the generation of vortex rings. This phenomenon 
is essential for large enough obstacles with the size about the healing length. In Bose-Einstein condensates (BECs) 
of rarefied gases the healing length can be relatively large and generation of vortices by small impurities becomes 
ineffective. In this case BEC remains superfluid for all velocities less than the minimal sound velocity corresponding 
to the long wavelength limit of the Bogoliubov dispersion law. For a supersonic motion of an impurity, the Cherenkov 
radiation of sound waves is the main mechanism of the appearance of friction and the corresponding "drag force" was 
calculated in 0, [1] . 

However, the detailed wave pattern generated by a moving impurity is also of a considerable interest. This problem 
became very topical in two-dimensional (2D) case in connection with the results of the experiment 0, 0| in which 
the waves were generated by the flow of a condensate expanding through an obstacle created by a laser beam. 
Since such an obstacle has the size comparable or greater than the healing length, the arising here wave pattern 
can be quite complicated. Already in numerical experiment Q modeling a similar situation it was noticed that the 
interference of sound (Bogoliubov) waves yields the wave pattern located outside the Mach cone. Analytic theory of 
such wave patterns was developed in @, [tj [13] • In many respects, this theory is analogous to the well-known Kelvin's 
theory of "ship waves" generated by a ship moving in still deep water, with the dispersion law for the surface water 
waves replaced by the Bogoliubov dispersion law for sound waves in BEC. Besides that, due to a large size of the 
laser beam, vortices or oblique dark solitons located inside the Mach cone can also be generated by a 2D flow of a 
BEC. The corresponding theory was developed in plElGl 

and was recently generalized [lij for a two-component 
condensate. However analogous theory for 3D flow has not been developed yet, although it is of considerable interest 
for understanding of the wave processes in BECs (see, e.g., [H). In this paper, we shall consider both analytically 
and numerically the 3D wave pattern created by a small impurity moving with supersonic velocity through a bulk 
Bose-Einstein condensate. 



II. STATIONARY WAVE PATTERN 



Dynamics of BEC of rarefied gases is described very well by the Gross-Pitaevskii (GP) equation 

ii/H + \M + (1 ~ M 2 )V> - Vi/> = , 



(1) 



FIG. 1: Coordinates defining the radius-vector r and the wave vector k. The latter one is normal to the wave crest line of the 
ship-wave which is shown schematically by a curve. 



which is written here in standard non-dimensional notation (see, e.g., 10]), so that the density of an undisturbed 
BEC with the repulsive interaction between atoms is equal to unity. We suppose that the external potential V is 
created by a point-like impurity moving with velocity U along x axis in negative direction, 

V(t) = V S(v + Ut). (2) 

The stationary wave pattern can be obtained for a supersonic velocity U which in our non-dimensional units with the 
sound velocity equal to unity means 

|U| = M > 1, (3) 

M being the Mach number. Assuming that the interaction with impurity is small, we can apply the perturbation 
theory [| S, G3 and linearize Eq. |T]) with respect to small disturbance (JvP of the wave function, = 1 + 5^. Then 
(5^ satisfies the equation 

iS^ t + ±A5f - {S^> + <5**) - V 6(r + Ut) = (4) 

which can be readily solved by the Fourier method and the resulting perturbation of density Sn — 6\^\ 2 = + 5^* 
is given by the expression (see, e.g., Eq. (20) in (Ioj |) 

f fc 2 e ikr d 3 k 

5n = Vo J (kU) 2 -/c 2 (l + /fc 2 /4)+ie(2^' ® 

where e = |e| • sgn(fc) is an infinitely small parameter, \s\ — > 0, which determines the rule of going around the poles of 
the integrand function. 

Since the wave pattern is axially symmetric with respect to the x axis, it is convenient to define the coordinate 
system so that the observation point lies in the (x, y) plane. Then vector r has the components 

r = (r cosx, rsinx, 0), (6) 

where x is the polar angle between r and x axis. Let the vector k lie in the plane making an angle <fi with the (x, y) 
plane. Then its components can be parameterized as 

k = (— fccosry, k sin 77 cos </>, k sin rj sin </>) , (7) 

where k cos r\ is the projection of vector k on the x axis. The geometrical meaning of the angles T],Xi an d H = 7r — rj — x 
is shown in Fig. 1. 

Substitution of Eqs. ([6]) and ([7]) into Eq. ([5]) and simple transformations cast this expression into the form 

V a f°° r k 2 smT 1 J (krsinr ] smx)e~ lkrcos,lcos x „ , 
6n = -*L L A^g-i dkd ^ w 
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where 

ko = 2y/M 2 cos 2 r]-l (9) 
and we have used the well-known integral representation 



J (z) = ±- e izcos *d<j> (10) 
27r Jo 

for the Bessel function. 

The integral in ([5]) can be estimated by the method similar to that used in Ref. First, reality of Sn enables one 
to represent Eq. (jSj) as a half-sum of this expression and its complex conjugate, and to make replacements k — > — k, 
e — > — e in one of integrals. As a result we get 

Vq r , f°° k 2 J a (krsmr]smx)e^ krcosrlcosx „ , , 

SU = 2^ J df]SmV P3fc2— *, (11) 

where fc-integration takes place over the whole k axis. The integrand function has two poles 

k = ±k + is' (12) 

located in the upper complex half-plane. Hence, we can calculate this integral, if we close the contour of integration 
by an infinitely large half-circle k = \k\e with either 0<9<norir<8< 2n, provided the contribution of these 
additional paths of integration vanish as |fc| — » oo. 

To analyze behavior of these integrals as |fc| — * oo, we use an asymptotic expression for the Bessel function, 



J o( z ) ~ y— cos^z- -J , 2>1. (13) 

Hence, Eq. (jTTJ) can be represented as a sum of two integrals with the integrand functions having the exponential 
factors 

exp[— ifcr cos(?7 ± x) — 7r/4]. (14) 

If cos(x ± rf) > 0, then this factor decays exponentially in the lower complex k half-plane, we close the contour by a 
lower half-circle, and since there are no poles inside this contour, the integral vanishes in this case. On the contrary, 
if 

cos(x ± r]) < 0, (15) 

then we close the contour of integration in the upper half-plane and both poles give non-zero contributions into the 
integral. Thus, we get 

Sn = — / k sin r\ ■ sin(fcr cos \ cos rj) ■ Jq (kr sin x sin if)d,Tj, (16) 
Jo 

where k is defined by Eq. (|9]), i.e. we have dropped out the index for convenience of notation. 

In the far-field region > 1 we replace the Bessel function by its asymptotic expression (fT3")) to obtain 



. Vo / k sin ?? 
on = \ / — ■ 



^ikr cos(x+7?)+7r/4 _^ ^ikr cos(x — if) — 7r/4 



dn. (17) 



7f V Sm X J0 

These integrals can be calculated by the method of stationary phase. The stationary point of the phase 

si = A:cos(x + rj) (18) 
is determined by the equation ds\/drj — which gives the relation between the angles x an d 

2M 2 

tan(x + r)) = — p - sin 2? 7 (19) 
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(l + fc 2 /2)tan X 

tanX= AP-(l + fcV2)- (20) 

This expression coincides with the results obtained in the case of 2D obstacle [tlQlj and satisfies the condition (|15p . 
On the contrary, the second term in Eq. lfl7|) with the phase S2 = k cos(x — rj) leads to the relation between \ an d V 
which is excluded by (JTSJ). Hence we take into account the first term only and reduce this integral into the Gaussian 
one around the vicinity of the stationary point. As a result we obtain the following distribution of the density in the 
wave pattern, 

Sn = ™ {[M^M* - 2) cos 2 „ + !][! + (4g/fl sh^]}^ 

7rr {[2M 2 cos 2 77-I] [l + (4M 2 /fe 2 ) cos 2?7 + (12Af 4 /fc 4 ) sin 2 2r/]}V2 1 VA " J ' v ; 

where x as a function of 77 is determined by Eq. (|20| and fc is defined by Eq. ([9]) . 

The geometric form of the wave crest surfaces can be easily found in the following way. Obviously, such a surface 
can be obtained by rotation of its cross section by the (x,y) plane around the x axis. Then we find from Eqs. ([6]), 
©, and (f2"0"l) the parametric formulae for the coordinates of this cross section: 

x = -g^ cos 77(1 — AT 2 cos 277) , y = ^ sin?7(2A/ 2 cos 2 77 - 1), (22) 

where s = kr cos(x + 77) is the phase constant along the crest line. These formulae are identical to ones obtained in 
2D case 0, [l(| which is natural since the Bogoliubov dispersion law for linear waves is the same for both two and 
three dimensions. However, the amplitude of waves as a function of the distance r and the polar angle \ (or 77) in 3D 
theory differs from that in the 2D case; now it decays with r as r -1 to satisfy the energy conservation law. 
As is clear from Eq. the wave pattern (|22[) corresponds to the range of the parameter 77 

- arccos(l/Af) < 77 < arccos(l/A/) (23) 

with the coordinates located inside the Mach cone defined by the relation 

1 , x 

sihxm = jj- (24) 

In particular, the small values of 77 correspond to the waves located in front of the obstacle, 

-~ ' . PM<-D V> vS (^L v . (26) 



2VM 2 - 1 4(A/ 2 - I) 3 / 2 ' ' w 2(Af 2 -l) 3 / 2 
i.e. the wave crest lines take here a parabolic form 



x(y) ~ ? + (M--lf^ 

(V) ~ 2VM^1 (2M 2 - l)s V • [ b) 



The boundary values 77 = ± arccos(l/M) correspond to the lines 



x 



±VM 2 - 1, (27) 

y 



i.e., far from the obstacle, they approach the straight lines parallel to the Mach cone (|24[) . In the region in front of 
the obstacle where y = z = 0, x < 0, we have 77 = 0, hence 



k = 2\jM 2 - 1 (28) 

and the wavelength 

A= ^ = _^ (29 ) 

k VAf 2 -l V ' 

is constant. Equation (|21[) reduces here to a simple formula 
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The formulae greatly simplify also in a highly supersonic limit and not too close to the Mach cone when M cos rj 3> 1. 
In this case Eq. (fT()|) yields 

X = 7T - 2,] -^AP (31) 

and even the leading order approximation \ — n ~ 2?/ gives good enough approximation in the most important region 
of the wave pattern. In particular, we get the expressions for the wave crest line 

x = — - (tan 2 n — 1), y — — tan 77, (32) 

2M v / j, y M h \ 1 

that is 

s M 9 .„„. 

x(y) = -7TTT + T^V (33) 



which is the limit M ^> 1 of Eq. ([26)) . The density disturbance (|2Tj) takes the form 

Vn 

=— cos M(r - x) , Mcosrj»l, (34) 
7rr 

and in front of the obstacle where x = — r = — |x| it corresponds to the limit M ^> 1 of Eq. (|3H)l . 

III. NUMERICAL SIMULATIONS AND DISCUSSION 

In our numerical simulations the GP equation 

Bib h 2 

ih-^ = -—Ai ) + V(r,t)^ + NU \i>\ 2 ^, (35) 



where 



Uq = 4irh 2 a s /m (36) 



is the effective interatomic coupling constant, a s being the s-wave scattering length of atoms, N is the number of 
atoms in the condensate, so that tp is normalized to unity, was transformed to non-dimensional units in the following 
way. We take some ao = ^Jh/muj x as a unit of length and w" 1 as a unit of time (if the BEC is confined in a parabolic 
trap then ao has a meaning of the "oscillator length" and uj x of the oscillator frequency along x axis) and introduce 

t = tuj x , r = r/a , i> = ■ a^ 2 , V = V/{mu 2 x al), g = 4Tra s N/a , (37) 
so that the non-dimensional GP equation takes the form 

i^^-^ + Viv^+g^ (38) 

with tildes omitted for convenience of the notation. 

In the current simulations the BEC was confined in a cubic box —10 < x,y, z < 10 and had practically uniform 
undisturbed distribution of density no — \ib\ 2 — 1.5714 • 10 -4 except for a narrow region at the boundary of the box. 
The other parameters have been chosen so that g = 8000, the sound velocity c s = y/gn^ — 1.1212, and the healing 
length £ = \/{-j2c s ) = 0.6307. The potential of the obstacle was represented by a spherical ball with the radius 
o-baii — 0.125 (which is less than the healing length) and the repulsive uniform potential equal to Vbaii = 100 inside 
the sphere. Velocity of the ball corresponds to the Mach number equal to M = 3. In our simulations we have used 
the method of lines with spatial discretization by Fourier pseudospectral method and time integration by adaptive 
Runge-Kutta method of order 2 and 3 (RK23). 

The resulting wave patterns is shown in Fig. 2. We have found that it is axially symmetric, as it was supposed, and 
the (x, y) cross sections of the wave crest surfaces agree very well with the analytical curves shown by dashed lines 
and corresponding to Eqs. I|22p. Oscillations of the density in front of the obstacle as a function of the x coordinate 
is shown in Fig. 3 and it is compared with the analytical expression (|30[) . Again good agreement is observed. Thus, 
the wave pattern located outside the Mach number is described quite satisfactory by the developed here theory. 



6 



-4 




X 



FIG. 2: (color online) Numerically simulated wave pattern generated by a spherical obstacle moving through a Bose-Einstein 
condensate. The parameters of the BEC are indicated in the text. The analytical wave crest lines are shown by dashed lines. 
They are plotted according to Eqs. (f22)l . 
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FIG. 3: (color online) Oscillatory structure in front of the obstacle obtained by numerical simulations (solid line) and analytically 
(Eq. (|30p ; red line). Position of the obstacle is shown by a blue dot. 



Since the size of the obstacle is much less than the healing length, there were no formation of vortex rings located 
inside the Mach cone. Just these structures attracted earlier much attention in the study of the loss of superfluidity 
in a subsonic motion of obstacles (see, e.g., [IB] and references therein) when stationary "ship waves" patterns do 
not exist — in subsonic case only time-dependent linear waves can be generated due to the switching on the obstacle 
potential [T7[ or a change of the obstacle velocity. Similar vortex-antivortex pairs are also generated in the 2D case 
where they align along straight lines as the velocity of the obstacle grows and above some critical value of velocity 
one can see the formation of oblique dark solitons attached at one their end to the obstacle and decaying into vortices 
at the other end. One may suppose that if the size of the obstacle exceeds the healing length, then in the 3D case 
formation of "conical solitons" would take place above some critical velocity. However, the study of this problem is 
outside the scope of the present paper. 

In conclusion, we have studied the formation of linear wave pattern generated by a 3D small obstacle moving with 
a supersonic velocity through a uniform condensate. Analytical formulae for the wave crest lines and dependence of 
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the amplitude of the density oscillations on the distance from the obstacle are confirmed by numerical simulations. 
This theory essentially extends previous calculations of the "drag force" and provides a more detailed picture of the 
process of Cherenkov radiation of Bogoliubov excitations in rarefied Bose condensates. 
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